function wint=interpolationphi_4agents(phi,jj, gammarfun, agentemp);
%%%trilinear interpolation%%%%%
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table elementsize
%jj=index(2);
lambda(1)=1/(1+gammarfun(1)+gammarfun(2)+gammarfun(3));
lambda(2)=gammarfun(1)/(1+gammarfun(1)+gammarfun(2)+gammarfun(3));
lambda(3)=gammarfun(2)/(1+gammarfun(1)+gammarfun(2)+gammarfun(3));
if lambda(1:3)>=lowerbound;
node=fix((lambda-lowerbound)/elementsize+1);
for g=1:3
   if node(g)==0
       node(g)=1;
   elseif node(g)==m
       node(g)=m-1;
   end
end
vv=ones(2,2,2,agentemp)*-9999;
xx=[node(1) node(1)+1]*elementsize;
yy=[node(2) node(2)+1]*elementsize;
zz=[node(3) node(3)+1]*elementsize;
for r1=0:1
    for r2=0:1
        for r3=0:1
            
            if node(1)+r1<=size(table,1) & node(2)+r2<=size(table,2) & node(3)+r3<=size(table,3)
            if table(node(1)+r1,node(2)+r2,node(3)+r3)~=0
            vv(r1+1,r2+1,r3+1,:)=phi(table(node(1)+r1,node(2)+r2,node(3)+r3),jj,:);
            end
            else
            vv(r1+1,r2+1,r3+1,1:agent)=-9999;    
            end
        end
    end
end


for g=1:agent
i1=vv(1,1,1,g)*(zz(2)-lambda(3))/(elementsize)+vv(1,1,2,g)*(lambda(3)-zz(1))/(elementsize);
i2=vv(1,2,1,g)*(zz(2)-lambda(3))/(elementsize)+vv(1,2,2,g)*(lambda(3)-zz(1))/(elementsize);
j1=vv(2,1,1,g)*(zz(2)-lambda(3))/(elementsize)+vv(2,1,2,g)*(lambda(3)-zz(1))/(elementsize);
j2=vv(2,2,1,g)*(zz(2)-lambda(3))/(elementsize)+vv(2,2,2,g)*(lambda(3)-zz(1))/(elementsize);

w1=i1*(yy(2)-lambda(2))/(elementsize)+i2*(lambda(2)-yy(1))/(elementsize);
w2=j1*(yy(2)-lambda(2))/(elementsize)+j2*(lambda(2)-yy(1))/(elementsize);
wint(g)=w1*(xx(2)-lambda(1))/(elementsize)+w2*(lambda(1)-xx(1))/(elementsize);
end
else
    wint(1:agent)=-9999;
end